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Abstract 

We study the decay of the probability for a non-Markovian stationary Gaus- 
sian walker not to cross the origin up to time t. This result is then used to 
evaluate the fraction of spins that do not flip up to time t in the zero tem- 
perature Monte-Carlo spin flip dynamics of the Ising model. Our results are 

compared to extensive numerical simulations. 
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Gaussian processes are amongst the most widely studied stochastic processes in various 
branches of science However, there are still simple but important questions associated 
with a Gaussian process that are nontrivial to compute. One such question is the following: 
consider a stationary Gaussian process X(t) with zero mean and a prescribed correlator. 
What is the probability P{(3) that X{r) does not cross the origin X = between r = 
and t = (31 This quantity, although simple and natural, turns out to be quite nontrivial to 
compute 0. Why is this so? A little thought shows that this quantity probes high order 
correlations in time of the dynamics and it depends on the whole history of time evolution 
of the system. In this Letter, we would like to address this question. It turns out that the 
solution to this problem has very wide applications in various other problems in Physics. For 
example, a related question arises naturally in the context of zero temperature Monte-Carlo 
dynamics in any spin system: what is the probability that a spin does not flip up to time 
t 0-|8[]. Similar questions also arise in the study of the fraction of lattice sites that remain 
unvisited up to time t by a random walker or by chemical species in a generic reaction 
diffusion system JTD|. 



In this Letter we restrict ourselves to one such application, namely the simple case of the 
Ising model. In the zero temperature dynamics of the Ising model, domains of opposite spins 
grow with time. At late times, the system is characterized by a single length scale (typical 
size of a growing domain) L(t) ~ t 1 ' 2 0. The fraction P(t) of spins that remain unflipped 



up to time t decays as P(t) ~ L(t)~ e for large time t ||, where 9 is a universal, dimension 
dependent, non-equilibrium exponent. Analytical computation of 9 seems to be extremely 
nontrivial. Even in d = 1, where the Glauber dynamics is exactly solvable, the calculation 
of 9 turned out to be a real tour de force, achieved recently by Derrida et al. ||. They found 
9\d = 3/4. However, their technique is special to d — 1 and seems impossible to extend to 
higher dimensions, where only numerical estimates of 9 are available @-§[. Therefore it is 
highly desirable to obtain an approximate analytical method to determine 9, that we now 
present. 

Let us start with a stationary Gaussian process X(r), with zero mean and a correlator 



(X(ti)X(t 2 )) = fiji — r 2 ). The probability P(/3) of not crossing the origin up to r = (3 
is expected to decay as exp(— 8(3) for large (3 (at least when /(r) decays exponentially for 
large r). We would like to calculate 8 since, as we will see later, 8 is related to the exponent 
8 of the Ising model. If the Gaussian process is Markovian, for which / is necessarily of 



8 = A exactly (see below also). For any other form of /, the process is non-Markovian (i.e., 
history dependent) and 8 is hard to compute. In fact, 9 depends very sensitively on the 
full function /(r) and not just on its asymptotic properties. Keeping the Ising problem in 
mind, we will restrict ourselves only to the class of non-Markovian processes for which (i) 
f'i^) 7^ an d (ii) /(t) ~ exp(— At) for large r. For convenience, we will then normalize /, 
setting /'(C^) = =Fl/2 after a proper rescaling of X, such that its Fourier transform satisfies 
u> 2 f(u) — > 1, for large u>. 

To illustrate the explicit history dependent nature of the non-Markovian process, it is 
useful to write its associated Langevin equation: 



where r](r) is a Gaussian white noise with (r](T)ri(r f )) = S(t—t'), and J is a causal (J(r) = 0, 
for t < 0) and integrable function. The history dependence is explicitly encoded in the kernel 
J. For J = 0, Eq. (1) describes a Markov process with /(r) = exp(— A|r|)/2A as stated 
above. In Fourier space, Eq. (1) amounts to X u = 77^(1 + J(uj))/(iuj + A), which allows us 
to relate the Fourier transform of / to that of J: f{oS) = |1 + J {uj)\ 2 / {uo 2 + A 2 ), with the 
correct large uo behavior, since J{oS) — * in this limit. 

We now proceed to a variational and perturbative calculation of 8, which will be tested 
by simulating Eq. (1), before applying these results to the spin flip problem. 

P(/3) can be written as the ratio of two path integrals, the first one Z\ over, say, positive 
trajectories X(r), the second one Z over unrestricted trajectories: 



the form /(r) = exp(— A|r|)/2A JTB|, it is possible to show by various methods PT5|— |T7]| that 




(1) 



2J x>0 VX(r)exp[-S} Z 1 
fVX(r)exp[-S] " Z Q 



(2) 
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where S = | Jq Jq X{t\)G{t\ — t 2 )X(t 2 ) drxdr%, and G(t\ — r 2 ) is the inverse matrix of 
f(r\ — r 2 ). 9 is then calculated from P{(3) by taking the limit, 9 = — Hindoo (3~ x lnP(/3). 
We impose periodic boundary conditions, X(0) = X(/3) for the paths, which should not 
affect the value of 9 in the limit of large f3. We notice that the Gaussian weight in Eq. 
(1) then becomes S = ^ J2 n G(uj n )\X(uj n )\ 2 , where G(u n ) = 1/ f(u n ) and u n = 2im/(3 are 
Matsubara frequencies. First consider a Markov process for which G(u) = u 2 + A 2 . We 
recognize the action in imaginary time (f3 is then the inverse temperature) of an harmonic 
oscillator of frequency A, S = fg C(X(r)) dr, with C(X) = \[{^f + X 2 X 2 }. Thus, P((3) is 
the ratio between the partition function of an oscillator with a infinite wall at X = and 
that of the same oscillator without the wall. For large (3, it goes as exp(—/3(Ei — Eq)), where 
Ei (Eq) is the ground state energy of the oscillator with (without) a wall. Thus, Eq = A/2, 
and Ei = 3 A/2, as the eigenstates of the problem with a hard wall at the origin are the odd 
states of the unrestricted oscillator. This gives the Markovian result 9 = A. 

For non-Markovian processes, unfortunately, S is no longer a classical action with an 
associated quantum problem. The denominator Zq can, however, still be computed exactly 
and we find after taking (3 — ► oo, 

As a check, one can verify that for an oscillator, for which G(u) = uj 2 + X 2 , one recovers E = 
A/2. The most difficult part is to evaluate the "ground state energy" Ei — — lim^oo \iaZi//3 
of the problem with a wall at the origin. One way of computing Zi will be to perturb 
around a classical action, for instance, that of harmonic oscillator of frequency ujq (or that 
of a particle in a box as we also did in fL7f). We adopt a variational method, by choosing a 
trial inverse correlator Gq(cj) = uj 2 + Uq, corresponding to that of an oscillator (with a hard 
wall at the origin) whose frequency Uq is going to be our variational parameter. We have 
the general variational inequality, Ex < 3uq/2 + lim^oo jj(S — Sq) w , where the average is 
performed using the action of the hard wall oscillator. The second term of the inequality 
requires evaluating the propagator (\X(u) n )\ 2 ) w of the hard wall oscillator, leading finally to, 
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(4) 



where the numbers c n can be evaluated using properties of Hermite polynomials and are 
given by 



(2n)! 



n\(2n - 1) 



(5) 



ir2 2n (2n + 1 

One can then differentiate Eq. (4) with respect to u to obtain an equation for uu Q , which 

(2) 

minimizes E\ . This equation can then be easily solved numerically. In principle, the value 
of Ei can be improved by summing higher terms of the cumulant expansion around the trial 
action ]17]. The superscript ^ denotes that we have kept only the first two terms of this 
expansion. To perform a systematic order by order cumulant expansion, one should also 
keep only the first two terms in E (even though E can be evaluated exactly to all orders 
from Eq. (3)). This gives 



E, 



(2) 



^0 



1 1 

— I / dx 
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G(xu ) 



1 



(6) 



u 2 (x 2 + 1) 

We can then define 9^ as 6^ = min a;o (£'i —E^f 1 ), remembering that 8 V = min WQ e[ 2 ^ —E , 
is an exact (presumably bad) bound of 9. 

When J(u) is small (see Eq. (1) and below), and using Eq. (4-6), one can perform a 
straightforward perturbative calculation around the Markov process J = to first order in 
K{uj) = J(u) + J(-u), 

2 1 



6 = A 



1 - -K(0) / K(x\)V(x) dx 

ix 2ix Jo 



(7) 



where V(x) = 4(x 2 + l)S(x) — 1, and S(x) is the same series that appears in Eq. (4). In 
fact, an infinite number of terms of the perturbation theory can be resummed by using a 
novel technique (for details, see JT7|), leading to, 



7T 



2tt 



6 = -JG(0) + — / W{G{uj)/lo z ) du 



n=l 



W(x) = J2 — ln(l + An 2 (x - 1)) - ln( 



x) 



(8) 
(9) 



This expression is valid provided G{uj)/uo 2 > 1, which is always the case for (and close to) 
a Markovian process. We have tested Eq. (8), by comparing its prediction to the direct 
simulation of Eq. (1), with /(r) = e exp(— r)/2 + (1 — e) exp(— 2r)/4, interpolating between 
two Markovian processes with A = 2 and A = 1. Note that e must be positive to ensure 
/ > 0, and that for e > 4/3, one can find u such that G{uj)/uj 2 < 1, and Eq. (8) is not 
reliable any more. The results for some representative values of e ^ 0, 1 (for which 8 = 2 
and 6 = 1) are displayed in Tab. 1, showing a very good agreement between Eq. (8) and 
numerical simulations. 

We now turn to the T = dynamics of the Ising model starting from a random (high 
temperature) initial configuration. We would like to show that calculating the fraction of 
unflipped spins up to time t in the Ising model reduces to calculating the survival proba- 
bility of a Gaussian process in the framework of Gaussian closure approximation (GCA), 



introduced by Mazenko fi3|J14]| . But before we make this connection, a few facts about the 
T = dynamics of the Ising model would be relevant. 

Following a quench to T = 0, interpenetrating domains of ±1 phases grow with time. A 
scaling theory has been developed to characterize the morphology of the growing domains 
According to this theory, at late times, the system is characterized solely by the linear 



length of a growing domain Lit). For the Ising model, Lit) ~ t 1 ^ 2 in all dimensions (however 
in d = 3 cubic lattice at T = 0, it seems that L(t) ~ t 1 / 3 due to lattice effects p2]Jl7|| , a 
fact which was underestimated in ||). Another prediction of the scaling theory relevant for 
us is that the on-site autocorrelation satisfies (S(t)S(t')) = F(L(t) /L(t')) for t > t' ^> 0, 



where F(x) ~ x for large x [^9],|TTJ- The exponent A is exactly 1 in d = 1, is close to 



1.25 in d = 2 (both from simulations |T^] and direct experiment [18[]) and close to 1.67 in 



d = 3 [|TTJ. GCA has been particularly successful in calculating this exponent as it predicts 
\id = 1, Am = 1.289 and X 3 d = 1.673, this approach becoming asymptotically exact for 
large d. Recently, we have extended GCA to the 9- state Potts model and calculated the 
g-dependent A PD| . 

Without entering into the details of GCA, we simply mention that this method assumes 



that the spin at position x and time t is essentially the sign of a continuous Gaussian 
stochastic variable m(x, t), which is physically interpreted as the distance to the nearest 
interface to the point x ||13|| . We thus see that in the framework of GCA, the probability 
that a spin does not change sign up to time t is equal to the probability that the associated 
Gaussian process m(t) does not cross the origin m = (from now we forget the label x 
as this probability does not depend on the considered site). The process m(t) at a given 
site is Gaussian whose correlator (m(t)m(t')) can be calculated using the GCA scheme. 
However, m(t) is not stationary since its correlation function depends explicitly on both 
t and t', and we need to use the following trick before using the results of Eq. (4-9). We 



define the variable X(t) = m(t)/y (m 2 (t)) which is also Gaussian. Its correlator turns out to 
satisfy (X(t)X(t')) = g(L(t)/L(t')), where once again g(x) ~ x~ x for large x. If we now set 
t = ln(L(t)), we are left to study the survival probability of a stationary Gaussian process, 
with correlator (X(t)X(t')) = f(r — r'), with /(r) = g(exp |r|). This completes the relation 
to the general problem studied in the first part of this Letter. For the Gaussian process, 
the survival probability P(r) ~ exp(— Or) ~ L{t)~ e . Thus we get 9 = 9 in the framework 
of GCA. Note that for an ordinary Brownian walker g(x) = x~ l l 2 , or /(r) = exp(— |r|/2), 
which leads to the well known result P(t) ~ t -1 / 2 . 

By means of GCA (see Jl7] as these calculations are not the subject of this Letter), we 
have computed the functions / for d = 1,2,3. For instance, in d — 1, we got the (exact) 



result, /(r) = J 2/(1 + exp(|2r|)). We see that / is not a pure exponential and is thus 
non-Markovian. Knowing /(r), we have computed the variational 6^ 2 \ and the associated 
frequency ujq. In Tab. 2, these results are compared to 9ns, extracted from numerical 
simulations of Eq. (1) (in Fourier space, then using inverse FFT), with the associated /. 
We also performed Ising simulations in d = 2 (800 x 800 lattice, 30 samples) confirming the 
results of @-§[), and in d = 3 (100 x 100 x 100 lattice, 60 samples). Note that the scaling 
of P(t) is found to be much better as a function of L(t) than as function of t pOjJ^. 



We note that the variational results are in good agreement with the simulations of Eq. 
(1), and with the Ising results in d — 1. However, the agreement in d = 2, 3 is not as good. To 



understand this, we have to remember that there are two possible sources of errors between 
the simulations of Eq. (1) and the Ising simulations: (i) the real variable X(t) ~ m(t) is 
not truly Gaussian and (ii) due to the use of GCA, the correlator / is not exact (except 
in d = 1, by accident). To illustrate the importance of the error due to (i), notice that 
although / is exact in d = 1, the Ising exact result 9 = 3/4 is bigger than the value 0.70 
obtained by simulating Eq. (1), which is perfectly reproduced by our theory. Note that, we 
got the exact bound for 6, 6 V = 0.736 in d = 1, whereas 9 V is a bad bound for d > 1. In 
d = 2, 3, the main cause of error is presumably (ii), as / given by GCA is found to decrease 
faster than in the Ising simulations for r < 2 |HJ. To check this, we computed 6 from Eq. 
(1) and 9^ using for / a fit of the correlator obtained in the Ising simulations, leading to a 
clear improvement (see Tab. 2). In d = 3, the numerical f(r) can be computed only up to 
a quite short time (r < 1.5), and 9 is quite dependent on the way f(r) is extended to large 
t. Indeed, for such short times, the known asymptotics /(r) ~ exp(— Ar) is not reached 
yet. Finally, one can easily understand why the theoretical estimates in d — 2, 3 seem less 
accurate than in d = 1 or in Tab. 1. It is possible to show [I7| that for the Ising model 



/ becomes less and less Markovian as d increases. An illustration of this is the fact that, 
as d increases, the relative difference between 9 and A strongly increases. In addition, the 
region for which G(uo)/uo 2 < 1 becomes wider as d increases, which is also a sign of strong 
"non-Markovianity" , and which prevented us from using Eq. (8) to evaluate 9. However, for 



large d, an alternate approach has been derived pl|JT7| which gives 9 in excellent agreement 
with numerical simulations. 

In summary, we have implemented variational and perturbative approaches to compute 
the survival probability of a Gaussian non-Markovian process, close enough to the Markovian 
limit associated to a quantum problem. These results were then used to calculate the 
fraction of unflipped spins in the zero temperature Monte-Carlo dynamics of the Ising model 
within the framework of GCA. Further details of these calculations, the application to the 
calculation of the g-dependent 9 for the (/-state Potts model p0|,^,|6| , numerics and other 



approaches will be presented elsewhere jLT]]. Finally the related question of the probability 



that the total magnetization never flips after a quench at T = T c |22| or at T = [17 is 
also the subject of further studies. 
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TABLES 



Tab. 1: 9 NS obtained from direct numerical simulations (NS) of Eq. (1) (T max ~ 1000, 
At ~ 10~ 3 , 100 samples, uncertainties are A9 = ±0.015) is compared to the result of Eq. 
(8) as a function of e, for /(r) = eexp(— r)/2 + (1 — e) exp(— 2r)/4. We have reported the 
point e = 4/3 (beyond which Eq. (8) must be regularized (see text)), which shows that 9 
can be less than the smallest inverse relaxation time (here 1). Note the remarkable accuracy 
of Eq. (8), up to its validity limit. 

Tab. 2: the four columns contain respectively, 1) 9i S i ng obtained in our Ising model numerical 
simulations for d = 2,3; 2) 9 NS obtained by simulating Eq. (1) (r maa , ~ 1000, At ~ 10 -3 , 
500 samples); 3) 9^ from Eq. (7-8); 4) the optimal frequency lu = 9^ which can be 
interpreted as the first term in the cumulant expansion. For d = 2,3, we also give the 
same results using a fit of / taken from the Ising simulations, instead of GCA. In d = 3, 
9 depends on the nature of the fit (see text) and we only give typically obtained values 
(estimated uncertainty of order 0.05). 
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TABLE 1 



£ 


0.10 


0.25 


0.5 


0.75 


1.33 




1.72 


1.47 


1.22 


1.09 


0.91 


^Theory 


1.74 


1.48 


1.24 


1.09 


0.91 



TABLE 2 





Rising 


8ns 


0(2) 




d= 1 


3/4 


0.70 ±0.01 


0.70 


0.61 


d = 2 (GCA) 


0.45 ± 0.01 


0.66 ±0.01 


0.58 


0.43 


d = 2 (fit of /) 




0.49 ± 0.02 


0.38 


0.27 


d = 3 (GCA) 


0.52 ±0.01 


0.70 ± 0.02 


0.47 


0.30 


d = 3 (fit of /) 




~ 0.55 


~ 0.4 


~ 0.2 
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